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TITLE 



TYPE 



NUMBER OF WORDS 
TEMPORARY STORAGE 
PARAMETERS 

STORAGE LOCATIONS 
3 



6 
7 



Solution of a System of Ordinary Differential. Equations 
(Originally Code 27) (SADOI Only) 
Closed - with one program parameter 



Loca t i on 
P 



p+1 



Program 
00 mF 
50 pF 



26 xF 



The first word of 
Routine F 1 - 11^ 
is at x« 



kl 

0, 1, 2, 3 

The locations 3 to 7 must contain the following parameters 
before and during the input of this subroutine* 
CONTENTS USE 

OOF 00 aF N(a + i) are the variables y. (i * 0, 
1, . .*, n - l) Originally the initial 
values are placed here* 

OOF 00 bF N(b + i) are the scaled derivatives, 2*hy! 
(=2 m hf 1 ), (i - 0, 1, 2, ,*., n - 1), 
calculated by the auxiliary subroutine* 
b > a + n-l 

Locations c+i, (i = 0, 1, «•*. n-l) 
are used as temporary storage for this 
subroutine. These locations must be 
cleared to zero before this integration 
subroutine is entered for the first time* 
c > b + n-l. 

n is the number of differential equations 
to be solved. 

d is the location of the first word of 
the auxiliary subroutine* 



OOF 00 cF 



OOF 00 nF 



OOF 00 dF 
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DURATION 



DESCRIPTION 



T = 7 + n(l5 + O.lm) + kt ms where 

T = time in milliseconds to perform one step of integration, 
t = time in milliseconds for the auxiliary subroutine. 
This subroutine will handle a set of n simultaneous first 
order ordinary differential equations, in which each deriv- 
ative is expressed explicitly in terms of the variables 



y 6 = f (y 0' y l' "•' y n-l } 
y l = f l (y 0' y l' •"' y n-l } 



^-1 = f n-l (y 0' y l' •"' Vl 5 

Any differential equation or set of differential equations 

to he solved must first be expressed in the above form before 

this subroutine can be applied. For example, the second order 

differential equation 

2 
y = w y 

must be written as two first order differential equations 



y =wy l 



T _ 



y[ = wy. 







where y, = y and y Q = y f / w - 



Each time this subroutine is called into use, it will carry 
out one integration step of length h. Each of the integrals 
y. (i = 0, 1, 2, ..., n-l) is replaced by its value at the end 
of a step of length h. In doing so, this subroutine employs 
an auxiliary closed subroutine which evaluates the functions 

f-, f . , f_, ..., f n from the given values of y . . The coder 

qj 2_' 2 n-l tf i 

must write this auxiliary subroutine for his individual 
problem since it defines the equations being solved and, this 
depends entirely on his specific problem. 

The purpose of the auxiliary subroutine is to calculate and 
store in locations b + i, (i = 0, 1, 2, ..., n-l), the 
quantities hf multiplied by a suitable scale factor 2 . 
h is the increment of the independent variable and,m is a 
positive integer to be chosen as large as possible without 
having any of the quantities 2 nf . exceed capacity anywhere 
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throughout the range of integration. The factor 2r is 
introduced to increase the accuracy of the integration 
subroutine- The variables, y. , must all "be scaled so 
that they are less than one throughout the range of 
integration before they are used in the auxiliary sub- 
routine. Also for maximum accuracy, one should store 
2 m h instead of Just h. This auxiliary subroutine must 
be located in a sequence of locations beginning with 
location d, where d is defined by the parameter S7« In 
integrating over one step, the integration subroutine will 
call in the auxiliary subroutine four times. 

This integration subroutine requires 3n arbitrary storage 
locations. The n consecutive locations a + i, (i = 0, 1, 2, 
..., n-1; a arbitrary), are used to store the variables y . 
It is in these locations that the initial values are to be 
placed. It is also in these locations that the final results 
are found. The n consecutive locations b + i, (i = 0, 1, 2, 
..., n-1; b > a + n-l), are used to store the scaled deriva- 
tives. 2 m h t ( = 2 m hf.), which are calculated by the auxiliary 
' y i 

subroutine. The n consecutive locations c + i (i = 0, 1, 2, 
..., n-1; c > b + n-1) are used for temporary storage by the 
integration subroutine. These locations will hold the 
quantities 2 m q. (See page 7). The numbers left in these 
locations at the end of an integration step are 3«2 times 
the roundoff errors of the quantities y . These numbers 
are taken into account during the following step and serve 
to prevent the rapid accumulation of roundoff errors . As 
a, result the effective numerical accuracy is m digits more 
than the capacity of the storage locations. Therefore, it 
is important that the locations c + i, (i = 0, 1, 2, ..., n-l) 
be cleared to zero hef ore the integration subroutine is 
entered. Otherwise, this integration subroutine will add 
spurious corrections to the variables . Thus before the 
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integration subroutine can be entered, the main routine 
must clear the temporary storage locations c + i to zero 
and set the initial values of the variables y in locations 
a + i. 
SUMMARY 

Supposing that, in the course of his routine, a coder 
has to solve a set of differential equations over a specified 
range given the initial value of the independent variable, 
a possible procedure would be the following: 

(1) Reduce the given set of differential equations to a 
set of n first order differential equations. 

(2) Calculate the initial values of the dependent 
variables, y . 

(3) Scale all the functions so that all the values y 
are less than one throughout the range of integration. 

(k) Choose a proper value of h (See note I). 

(5) Choose m properly. v 

(6) Determine the parameters to be placed in S3 - S7, 
observing that a < b < c. 

(7) Write an auxiliary subroutine which evaluates the 
functions 2 nf and stores them in locations b + i. 

(8) Make certain that the main routine sets the scaled 
initial values in locations a + 1, and clears the temporary 
storage locations c + i to zero before the integration 
subroutine is entered. 

With respect to the solution of a set of differential 
equations, a program can be broken up into four parts: 

(1) ■* Locations 3 through 7 which contain the parameters, 

(2) The main routine, 

(3) The integration subroutine (Code F 1 - 11*0 
(k) The auxiliary subroutine. 
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THE IHDEPEHDEflT VARIABLE 

If the independent variable x occurs in the functions 
f or if it is required during an integration as an index, 
then it must be obtained by integrating the equation 
x T = 1. The independent variable x is then treated as an 
additional dependent variable, for which the auxiliary- 
subroutine has to provide the quantity 2 m hx t = 2 m h. However, 
this latter quantity may be planted at the beginning of the 
integration in the appropriate location (e.g. in location b) 
and left there, so that the auxiliary subroutine is relieved 
of the task. If the independent variable does not appear 
in any of the f ' s but is merely wanted for indication 
purposes, it is quicker to use a simple counter in the main 
routine . 

NOTES 

I) Accuracy: The truncation error in one step is of the 
5 
order of h . Ordinarily, that is for a small set of well 

behaved equations, its magnitude is about 10 h ; for 

large sets or difficult equations it may be greater. Over 

the range of integration this error will amount to about 

h /100. Roundoff errors accumulate at a rate corresponding 

to the keeping of (39 + m) binary digits . The choice of 

the length of the increment h is governed largely by the 

accuracy desired. An increase in the length of h will 

result in a decrease in accuracy and in operating time. 

Likewise, a decrease in the length of h will result in an 

increase in both accuracy and operating time. However, 

no further increase in accuracy can be gained by choosing 

# h < 2 because of the introduced truncation error. But, 

if the functions are very sensitive to variations in y , 

or if the number of equations is very large, smaller steps 

will probably be necessary with, of course, a corresponding 

increase in the time required. How, the process used in the 
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integration subroutine is a fourth order one. Thus, 
l/l5 of the following difference, 

(the value of y calculated using an interval 

of length h) 
-(the value of y calculated usin^g an interval of 

length 2h) 
is an approximation of the error. 

II) Adjustment of the increment $$ There exist 
essentially two ways of adjusting the increment 

a) One may double or halve the inc%ment "by varying 
the value of m in the main routine. ^Th^ ^V De done over 
the complete range of integration or j^us^ over part of it. 
When only the parameter 00 mF in the ifyiSt "between the main 
routine and the auxiliary subroutine is el^anged to 00 (m+l)F 
and the auxiliary subroutine is unalter^i the length of the 
increment is halved. Likewise, when onl$| the parameter 00 mF 
is changed to 00 (m-l)F the length of the. increment is 
doubled, tfhe auxiliary subroutine is not Altered since 
2%. = 2 m+ Ji/2t If one adjusts the increment over the complete 
range, adjustijig only the value of m is sufficient. However, 
if one wishes fco adjust the length of the increment within 
the range of integration, one must also adjust all the 
quantities in locations c + i. Otherwise, one will introduce 
roundoff errors in y. of the magnitude, 

1 -ko 

2 (old value of h - new value of h) x 2 
Now by also doubling the quantities in c + i when one halves 
the increment one will introduce no roundoff error. Similarly, 
by halving the quantities in c + i when one doubles the 
length of the increment, one will introduce no roundoff 

error. If one clears the locations c + i, one introduces 

-1*0 
roundoff errors of magnitude 2 
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b) One may alter the length of the increment in any 
ratio by adjusting the scaling factor 2% in the auxiliary 
subroutine. Here also one may adjust the length of the 
increment within the range of integration. Now it is not 
necessary to adjust the quantities in c + i. If, however., 
2 li becomes small, then roundoff errors are introduced by 
inaccuracies in the auxiliary subroutine. Thus one should 
not keep 2% small when integrating over large ranges unless 
the loss of accuracy and time does not matter. 

Ill) Often it is desired to evaluate functions involving 

expressions like sin x or J (x). These expressions can be 

m 

evaluated by solving extra distinct differential equations 

along with the desired ones. For example, 

2 2 
d /dx (sin x)/2 = -(sin x)/2. 

Thus we can evaluate (sin x)/2 by using the extra pair, of 
equations 

y JLi = 2,D hy y' = -2 1B hy 

'n+1 v n •'n ^n+1 

and suitable initial conditions . 

METHOD USED FOR INTEGRATION IN THE ROUTINE 
Given a set of differential equations, 

y i = f i ( y o' y l' y 2 , ""* y n-l'' ^ = °* ls 2 > '•> n - 1 ) 

The process used in the integration is defined by the 
following equations 

k ij =2lDhf i froj' y lj> — y n -lj> 

P 1,J + 1 * <Vl + ^ ^j-Vi^ 

y i,j+l = y i,j + 2 r i,j + i 

a, . , = q . . + 3r + (C . - 1) k, . _ 
i*J+l ij i,j x j i*J+l 

with the following table of values 
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J 


J+l 


B. 
J 


C. 





-1/2 


2 


1/2 


1 


"(1/2) 1/2 


1 


(1/2 ) 1/2 


2 


(1/2 ) 1/2 


1 


-d/2) 1 / 2 


3 


-5/6 


2 


1/2 



Of the double subscripts ifsed in the above equations, the 
first subscript, i, indicates which variable is being 
considered, and the second Subscript, j, indicate which of 
the four parts of one step it being performed. The auxiliary 
subroutine evaluates the quantities k„ . . In the above 
equations, only the quantities q. , and y. . are carried 

1,4 1,4 

over from step to step. The quantities r. are calculated 
in the course of one step; they are not carried directly 
from step to step. When j = 4, we replace it by zero, 
increase i by 1, and terminate the step. 

For one step, the sequence of operations is as follows: 



J = o 


i = 0, 1, 2, ..., n-1 


j = 1 


i = 0, 1, 2, ..., n-1 


J = 2 


i = 0, 1, 2, . . ., n-1 


J = 3 


i = 0, 1, 2, ..., n-1 
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LOCATION ORDER 


__ _. BOTES 


PAfiS 1 








00 K(P1) 














00 S3 














00 g4 




IK 


| 






00 S5 












00 86 














01 29K j 













L5 5F 












L4 6f 












1 


00 20F 
46 571F 












2 


L5 4f 














LO 3F 


1 








3 


42 573F 
00 20F 












V 


46 573F 
L5 5F 












5 


LO 4f 
42 574F 












6 


00 20F 
46 574f 












7 


26 93F 
00 F 












8 


64 f 
00 33L 




c+n 








9 


80 S3 

00 S3 




a 

a 








■10 


00 F 




b-a 






i 

! 
i 


11 
12 


00. F 

00 F 1 

00 F 

26 1469N 

01 K 




b-a 

c-b 






; 




l 


S5 F 
46 8L 
46 HL 




Set shift addressee 
= m 






' 


j L4 3L 


1 




J 





F 1 



LOCATION 


ORDER 


NOTES PAGE 2 


2 


42 22L 
22 21L 




_| Set link address 


3 


L5 (q^)* 

40 IF 


From 21 
By 20 




4 


L5 (ic^jF 

40 F 


By 18 




5 


LO (1)F 


By 23 




6 


LO IF 
40 3F 




Ki- B i^ 




50 (Aj)F 


By 24 




7 


7J 3F 






8 


L4 3F 
10 (m)F 
40 3F 


By 


^ij-Vu 5 V 1)2 " 


9 


L4 ( yiJ )F 


By 17 






10 


40 (y. . . )F 

L5 3F 


By 17 




Step y 


11 


50 2F 
00 (m)F 
40 3F 


By 1 


Form r. . _ 
i*J+l 


12 


50 (Cj)F 


By 25 




13 


7 J F 
LO F 
L4 IF 




C. k. . - k. . 


14 


L4 3F 






15 


L4 3F 
L4 3F 




q ij + 3 r i,j+l 


16 


k0 kl,J + l> F 
L5 9L 

L4 13L 


By 19 




17 


42 9L 
46 9L 






Increase all addresses depending 


18 


L4 39L 
46 4l 






on i by 1 



Fl 



LOCATION 


ORDER 


NOTES PAGE 3 Fl 


19 


l4 4ol 

42 15L 






■ ■ 1 


20 


46 3L 
LO 37L 






until i = n 


21 


36 3L 










L5 (33 )L 


By 22, 






22 


k2 21L 


From 2 f 




Increase j from to 3 and then leave 




32 ( )F 


By 2 




by link 


23 


46 5L 
10 10F 






Adjust addresses which depend on j 


24 


L4 18L 
42 6L 








25 


46 12L 






■ 




50 25L 




Call in auxiliary subroutine 


26 


26 S7 








41 2F 


• 


Clear 2F so that it can be used as zero. 


27 


L5 38L 




. 




26 1?L 




Start new i cycle « 


28 


40 F 






i 


00 F 


1/2 


C 0' C 3 


29 


NO F 








00 F 


-1/2 


A 


30 


40 F 






! 


oo 2071 0678 1186 J 


1 jf2 C 1 , k 2 


i 


80 F 






i 


00 2928 9321 88l4 J 


-1 //¥ k ±s c 2 


32 


80 F 








00 1666 6666 666j j 


_-5/6 A^ 


33 


LJ 1025F 


-11, 1 








06 1058L 


25,34L 




Expressed in units of 2~ 9 , 2*" 19 , 2~ 29 , 2* 59 . 


34 


LJ 3074F 


-9 2 




Addresses used to set addresses to refer to 




06 3107L 


27, 35L 




the constants A„, C. and make the address in 


35 


LF 2F 


-8, 2 




5L, 1 or 2 according as B. 2 or 1, and to 

J 
stop the address in 21, dependent on j, and 




06 2084L 


26, 36L 




36 


LJ 1Q25F 


-11, 1 




to stop when positive. 




07 37L 


28, 37L 








0141K 







